Europhysics Letters 



PREPRINT 



Free-Energy Barriers in the Sherrington-Kirkpatrick Model 

Elmar Bittner and Wolfhard Janke 

Institut fur Theoretische Physik and Centre for Theoretical Sciences ( NTZ). 
Universitat Leipzig, Augustusplatz 10/11, D-04109 Leipzig, Germany 



PACS. 02.70.Uu - Applications of Monte Carlo methods. 
PACS. 75.10.Nr - Spin-glass and other random models. 



Abstract. - The Sherrington-Kirkpatrick spin-glass model is investigated by means of Monte 
Carlo simulations employing a combination of the multi-overlap algorithm with parallel tem- 
pering methods. We investigate the finite-size scaling behaviour of the free-energy barriers 
which are visible in the probability density of the Parisi overlap parameter. Assuming that the 
mean barrier height diverges with the number of spins N as N a , our data show good agreement 
with the theoretical value a = 1/3. 



Despite three decades of research the nature of the "glassy" low-temperature phase of 
finite-dimensional spin-glass systems remains a major open problem in statistical physics. It 
is still unresolved whether the replica symmetry-breaking theory or the phenomcnological 
droplet picture yields the correct description (for reviews, see refs. [1-4]). Even at the mean- 
field level, only very recently a mathematical proof [5] of Parisi's replica solution [6] for the 
Sherrington-Kirkpatrick (SK) model [7] was given. 

In the thermodynamic limit the frozen phase of the mean field spin glass shows many stable 
and metastable states. Such a feature is the consequence of the disorder and the frustration 
characterizing spin glasses in general and leading to a rugged free-energy landscape with 
probable regions (low free energy) separated by rare-event states (high free energy). But 
also for finite systems the free-energy landscape shows an intricate, corrugated structure. 
Therefore, it is hard to measure the free-energy barriers by means of conventional Monte 
Carlo simulations directly. The aim of this letter is to study the free-energy barriers of the SK 
mean field spin-glass model using the multi-overlap Monte Carlo algorithm [8] and to compare 
our results with previous findings for finite-range spin glasses [9]. For analyzing the barriers, 
we used precisely the same method as introduced in ref. [9] for the Edwards- Anderson (EA) 
nearest-neighbor model [10], where the scaling of the mean barrier height with the number of 
spins was found to clearly deviate from the theoretical mean-field prediction in both three and 
four dimensions. To exclude the possibility that this deviation could, in principle, be caused 
by employing different definitions in the theoretical and computational work, it is important 
to apply precisely the same numerical procedure to the SK model. 

The Hamiltonian of the SK model reads 
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where Sj = ±1, i = 1,. . . , N, and the are independent random variables with a Gaussian 
distribution of zero mean and variance A^ -1 , N being the numbers of spins. The critical 
temperature of the infinite system is T c = 1. 

The fact that there is no parametrization of the relevant configurations by a conventional 
thermodynamic variable led us to use the Parisi overlap parameter [6] , 



1 N 

q= /vE 



v 2-i (1 M a) , (2) 

1=1 

as an order parameter, where the spin superscripts label two independent (real) replicas for 
the same realization of randomly chosen exchange coupling constants J = Jij. For given J 
the probability density of q is denoted by Pj(q), and the function P(q) is obtained as 

^) = [5(?)]av = ^E F ^)' ( 3 ) 
^ J J 

where # J is the number of realizations considered. For a given realization of J the nontrivial 
(i.e., away from q = ±1) minima are related to the free-energy barriers of this disordered 
system J. We are, therefore, interested in the whole range of the probability density Pj{q). 
Conventional, canonical Monte-Carlo simulations are not suited for such systems because 
the likelihood to generate the corresponding rare-event configurations in the Gibbs canonical 
ensemble is very small. This is overcome by non-Boltzmann sampling [11, 12] with the multi- 
overlap weight [8] 



wj(q) = exp 



i<j 



(4) 



where the two replicas are coupled by Sj in such a way that a broad multi-overlap histogram 
Pj Uq (q) over the entire accessible range — 1 < q < 1 is obtained. When simulating with the 
multi-overlap weight, canonical expectation values of any quantity O can be reconstructed by 
reweighting, 

(O)T = (O eMSj))j/(eMSj))j ■ (5) 
Ideally the weight function Wj = exp(Sj) should satisfy 

Pj uq (q) = Pj n Wj = const. , (6) 

i.e., it should give rise to a completely flat multi-overlap probability density Pj Uq (q). Of 
course, Pj(q) is a priori unknown and one has to proceed by iteration. An efficient way to 
construct the weight function Wj is to use an accumulative recursion, in which the new weight 
factor is computed from all available data accumulated so far, for details see refs. [13] and 
[14]. The multi-overlap algorithm combined with this recursion allows an almost automatic 
simulation of the SK model. 

Let us now turn to the description of the Monte Carlo update procedure used by us. 
We combined the multi-overlap algorithm [8] as described above with the parallel tempering 
optimized Monte Carlo procedure [15] to overcome as many as possible hidden barriers in the 
rugged phase space. We studied systems with N = 32, 64, 128, 256, 512 and 1024 spins. For the 
parallel tempering procedure we chose a set of 32 temperature values in the range T = 1/3—1.6 
for all of our systems apart from the largest, where we used 64 temperature values for the 
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same temperature interval. Once the entire range of q for all temperature values was covered, 
the accumulative recursion for the weight functions was stopped. Due to large differences in 
the free-energy landscape for different disorder realizations Sf, the number of recursion steps 
varied for different J . After the weight functions were constructed, they were kept fixed and 
we took about 100 000 measurements, with five sweeps between the measurements. A sweep 
consisted of N spin flips with the multi-overlap algorithm and one parallel tempering update. 
We recorded time series of the overlap parameter q for five different temperature values and 
the canonical Pj{q) distribution for all temperature values. To average over the disorder we 
used 1000 realizations of the disorder for N < 512 and 100 for N = 1024. 

For each of these samples J we computed the barrier autocorrelation time tb by employing 
the same method as Berg et al. used for the ±J EA Ising spin-glass model [9]. For clarity, 
we recall the basic idea here. The free-energy barrier Fg for a given Pj(q) is defined through 
the autocorrelation time of a one-dimensional Markov process which has the canonical Pj(q) 
distribution as equilibrium state. The transition probabilities Tjj are given by 



T 



1 - io 2 ,i 
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W 3 ,2 









(7) 



where Wij (i ^ j) is a probability a la Metropolis to jump from state q = qj to q = qi 
{ qi = i/N,ie[-N,-N + 2,...,+N}), 



1 ■ f-. Pj{Qi 
— mm II,- 



(8) 



The transition matrix T fulfils the detailed balance condition (with Pj), and as a consequence 
it has only real eigenvalues. The largest eigenvalue (equal to one) is non-degenerate, and the 
second largest eigenvalue Ai determines the autocorrelation time of the Markov chain, 



1 



T B,J 



iV(l-Ai) 

The associated free-energy barrier for realization J is defined as 



B,J 



In (tb 



J) 



(9) 



(10) 



Note that the definition of the autocorrelation time (© takes only barriers in q into account, 
but not other barriers which may well exist in the multidimensional configuration space. 

For each temperature value we performed least-squares fits of the finite-size scaling (FSS) 
ansatz Fb = [Fb,j]d, v — cN a which corresponds to the exponential FSS behaviour 



tb oc e 



cN a 



(11) 



The results from these fits are depicted in fig. ^ and are consistent with previous results in 
the literature [16-23] using numerical and analytical methods. The horizontal line in fig. ^ 
indicates the theoretical value a = 1/3 of ref. [21]. The three data sets show fits with different 
lower bounds of the fit range iV m i n , while the upper bound was always our largest system 
N = 1024. From these fits we sec a strong finite-size effect for T — > T c = 1. At lower 
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Fig. 1 - Dependence of the exponent a in the ansatz Fb = cN a on the lower bound of the fit range 
[iVmin, 1024] as a function of temperature. The horizontal line indicates the theoretical value a = 1/3. 



temperatures we find a linearly increasing deviation from the theoretical value. This is also a 
finite-size effect, because the slope of the deviation becomes natter when increasing the lower 
bound of the fit range and there is no physical reason for a change of behaviour of the barrier 
autocorrelation time in the glassy phase. 

One possible explanation for this deviation from the theoretical value is the lack of self- 
averaging of the finite volume Parisi overlap parameter distribution Pj in the SK model [24]. 
To check the non-self-averaging of our numerical data, we analyzed our free-energy barriers 
relying on the (empirical) cumulative distribution function D(x), which is defined for a set of 
sorted data, e.g., the free-energy barriers Fb,i < Fb,2 < ... < Fb, u , by 

- - =- < D{x) <- + ^~ for Fb i <x< F B , l+l , (12) 
n 2n n 2n 

where we use a straight-line interpolation in between. A nice way to test the non-self-averaging 
is to look at the peaked distribution function [9] : 

(D(x) for D(x) < 0.5 , 
Uq{X > ~ \l- D(x) for D(x)>0.5 . [ 6) 

This function has a peak at the median x me d of the data with Dq = 0.5. For self-averaging 
data x, the function Dq collapses in infinite volume to Dq(x) = 0.5 for x = x and otherwise, 
where x is the mean value. For non-averaging quantities the width of Dq stays finite. In 
fig. [21 we show the peaked distribution function in units of the median value which also scales 
approximately with TV 1 / 3 . We observe a very weak finite-size dependence and therefore clearly 
a non-self-averaging behaviour of the free-energy barriers. This result is in agreement with 
theoretical calculations [24]. The width of Dq increases with decreasing temperature, i.e. the 
sample to sample variations become more pronounced for lower temperatures. 

We already mentioned that the distribution of the free-energy barriers becomes broader 
for low temperatures. Now let us have a closer look at the distribution itself, therefore we 
only analyze the system sizes where we have 1000 different disorder realizations. In recent 
work, Dayal et al. [25] have found that the tunnelling times of their flat-histogram sampling 
simulations of the 2D ±J EA Ising spin glass are distributed according to the Frechet ex- 
tremal value distribution for fat-tailed distributions. In general, extreme-value statistics can 
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Fig. 2 - Probability distribution function Dq for the free-energy barriers at T — 1/3 in units of their 
median value. 

be classified into different universality classes [26,27], depending on whether the tails of the 
original distribution are fat tailed (algebraic), exponential, or thin tailed (decaying faster then 
exponential). Assuming that the free-energy barriers of the SK model are, as well as the tun- 
nelling times of the EA spin glass, distributed according to an extreme-value distribution, we 
use the integrated probability density of the generalized extreme value distribution (GEV), 



Ft,^( x ) = ex P 



-i/C 



(14) 



for 1 + £(x — fi)/cr > 0, to fit our data. Here, £ is the shape parameter, and [i and a are 
related to the A-dependent mean and variance of the distribution, respectively. We find that 
the free-energy barriers show fat tails for T < T c with shape parameter £ > 0, i.e., a Frechet 
distribution. In fig. we plot the resulting fits and data of the probability density for different 
temperatures below the glass transition and find that the tails become fatter and fatter as 
the temperature goes to zero. The histograms for low temperatures show deviations from the 




Fig. 3 - Density of free-energy barriers Fb for TV = 256 at different temperatures. The inset shows 
the densities for T = 1/3 for different numbers of spins. 
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Fig. 4 - Temperature dependence of the parameters /i and a of the Frechet distribution for N = 512. 
The inset shows the size dependence of fx and a for T = 0.394. 



Frechet distribution for small values of Fb, so a much larger number of disorder realizations 
would be needed to determine both tails of the distribution properly. We determined the 
parameters fi, a and £ for different temperatures and found that \x grows logarithmically and 
a linearly with inverse temperature X/T, whereas £ stays more or less constant at £ ~ 0.33. 
As an example we show in fig. 0] the results for N = 512. If we keep the temperature fixed 
and look at the size dependence of the distribution, we find that for a larger number of spins 
the distribution becomes broader, c.f. the inset of fig. To quantify this behaviour we use 
the scaling relations \i oc N a » and a oc N a " , which lead to a M w 0.31 and a a w 0.25 for 
our lowest temperatures, see the inset of fig. 0| The exponent of the location parameter fi is 
closely related to the scaling exponent of the mean barrier height, therefore should be 1/3, 
which is in good agreement with our result. We find a weak temperature dependence of the 
exponents a M and a„ with positive and negative slope for increasing T, respectively. 

As a short remark we want to mention that we also looked at the correlation between the n 
largest eigenvalues A„ of the interaction matrix Jy, which is a real symmetric matrix, and the 
largest free-energy barrier Fb for this disorder realisation. And we find for T < T c a very weak 
but non- vanishing correlation between the largest eigenvalue Ai and Fb- Unfortunately, this 
correlation becomes weaker and weaker for larger systems and also for smaller temperatures. 
All other correlations vanish much faster, such that these correlations cannot be used to 
predict the magnitude of the free-energy barrier based only on Jjj . 

To conclude, we found that the free-energy barriers of the SK model are non-self-averaging 
and distributed according to the Frechet extremal value distribution. These particular features 
were also found for the EA nearest-neighbor model and such similarities support the position 
that the Parisi replica symmetry breaking solution of the SK model is the limit of the short- 
range model on a lattice in dimension d when d — > oo, with a proper rescaling of the strength 
of the Hamiltonian. On the other hand, we also found that the free-energy barriers diverge 
with the theoretically predicted value a = 1/3, which is in contrast to the results for the 
EA model in three and four dimensions [9]. Of course, one reason for this discrepancy could 
possibly be that finite-size effects, caused by the relative large temperature used in ref. [9], 
are too strong or the lattice sizes are too small to see the real asymptotic behaviour. Still, our 
results also support the complementary position, that the Parisi replica symmetry breaking 
solution sheds little light on the thermodynamic structure of the EA model [28]. With these 
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oppositional results we are not able to rule out one of the two scenarios. But we can state that 
the method to determine the free-energy barriers proposed by Berg et al. [9] leads to correct 
results for the SK model and is very useful for systems with many barriers in the free energy. 
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